FUNCTION comovingdistance_integrand,z

COMMON comovingdistance_internal,omeg0,lam0

return,1./sqrt((1.+z)^2*(1.+omeg0*z)-z*(2.+z)*lam0)

END 

;**************************************************************************


FUNCTION comovingdistance,redshift,omega0,lambda0,hubble0,help=help

COMMON comovingdistance_internal,omeg0,lam0

; Returns comoving distance (Mpc) to redshift z.
; Assumes only matter and vacuum energy important.
;
; Written by Joop Schaye

IF keyword_set(help) THEN BEGIN
    print,'SYNTAX: result = rcomov(redshift,omega0,lambda0,hubble0)'
    print,'Default: returns comoving distance (Mpc) to redshift z'
    return,0.
ENDIF

IF lambda0 EQ 0. THEN BEGIN
    comdist = omega0*redshift + (omega0-2.0)*(sqrt(1.+omega0*redshift)-1.0)
    comdist = comdist/omega0^2/(1.+redshift)
    comdist = comdist * 2.0 * 2997.9/hubble0
ENDIF ELSE BEGIN ; Use Carroll & Press 1992, eq 25.
    omeg0 = omega0
    lam0 = lambda0
    integral = qsimp('comovingdistance_integrand',0.,redshift)
    omegak=1.-omega0-lambda0
    IF omegak EQ 0. THEN BEGIN
        comdist = integral
    ENDIF ELSE BEGIN
        IF omegak GT 0. THEN comdist = sinh(sqrt(abs(omegak))*integral) $
        ELSE comdist = sin(sqrt(abs(omegak))*integral)
        comdist = comdist / sqrt(abs(omegak))
    ENDELSE 
    comdist = comdist * 2997.9/hubble0
ENDELSE 

return,comdist

END 
